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We  present  a  methodology  for  the  efficient  calculation  of  the  shock  Hugoniot  using  standard 
molecular  simulation  techniques.  The  method  is  an  extension  of  an  equation  of  state 
methodology  proposed  by  Erpenbeck  [1992,  Phys.  Rev.  A,  46,  6406]  and  is  considered  as  an 
alternative  to  other  methods  that  generate  Hugoniot  properties.  We  illustrate  the  methodology 
for  shocked  liquid  N2  using  two  different  simulation  methods:  (a)  the  reactive  Monte  Carlo 
method  for  a  reactive  system;  and  (b)  the  molecular  dynamics  method  for  a  non-reactive 
system.  The  method  is  shown  to  be  accurate,  stable  and  generally  independent  of  the 
algorithm  parameters.  We  find  excellent  agreement  with  results  calculated  by  other  previous 
simulation  studies.  The  results  show  that  the  methodology  provides  a  simulation  tool  capable 
of  determining  points  on  the  shock  Hugoniot  from  a  single  simulation  in  an  efficient, 
straightforward  manner.  Further  applications  and  extensions  of  the  method  are  briefly 
discussed. 


1.  Introduction 

The  behaviour  of  materials  under  conditions  of 
extreme  temperature  and  pressure  is  of  significant  interest 
in  many  fields  of  physics  and  fluid  science  [1-4].  Of  special 
interest  are  energetic  materials,  a  class  of  materials 
of  critical  industrial  and  military  importance.  These 
materials  exhibit  chemically  and  physically  interesting 
behaviour  when  exposed  to  extreme  temperatures  and 
pressures.  In  particular,  when  subjected  to  shock, 
energetic  materials  often  undergo  rapid  reactions  that 
produce  a  heterogeneous  mixture  of  chemical  species 
that  are  accompanied  by  huge  energy  releases  and  can 
produce  pressures  up  to  several  hundred  GPa  and 
temperatures  exceeding  10  000  K.  For  a  sufficiently 
strong  shock,  a  supersonic,  self-propagating  reaction 
wave  known  as  a  detonation  can  be  initiated.  Unfortu¬ 
nately,  the  extreme  conditions  along  with  the  short  time 
and  length  scales  over  which  a  detonation  occurs  poses 
considerable  experimental  challenges  in  characterizing 
the  material  behind  the  detonation  front.  Therefore  a 
concerted  effort,  which  combines  experimental,  theo¬ 
retical  and  simulation  approaches,  is  essential  for 
furthering  our  understanding  of  such  shocked  systems. 
Advances  in  experimental  capabilities  provide  us  with 
crucial  property  data,  while  the  continuing  development 
of  accurate  equations  of  state  has  allowed  reasonable 
predictions  of  various  shock  properties  [5,  6].  Similarly, 
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the  development  of  novel  methods  to  simulate  these 
complex  systems  has  been  the  focus  of  research  efforts 
and  has  recently  led  to  the  invention  of  some  uniquely 
effective  simulation  tools  [7^4].  These  classical  simula¬ 
tion  methods  can  be  implemented  irrespective  of  rate 
limitations,  the  production  of  huge  energy  releases,  or 
extreme  thermodynamic  conditions. 

The  Hugoniot  curve,  a  commonly  calculated  property 
in  shock  and  detonation  science,  reveals  many  proper¬ 
ties  of  shocked  materials  and  knowledge  of  which  is 
critical  to  the  design  of  new  materials  and  application 
platforms.  This  curve  consists  of  the  set  of  (PVT)  points 
for  which  the  Hugoniot  expression, 

Hg  =  E-E0-\(P  +  P0)(V0-V),  (1) 

is  zero.  In  equation  (1),  E  is  the  specific  internal  energy, 
P  is  the  pressure  and  V  —  1/p  is  the  specific  volume 
(p  is  the  specific  density).  The  term  specific  refers  to  the 
quantity  per  unit  mass,  while  the  subscript  ‘o’  refers  to 
the  quantity  in  the  initial  unshocked  state. 

Presently,  three  approaches  exist  for  calculating  the 
shock  Hugoniot  states  from  classical  molecular  simula¬ 
tion,  each  with  their  own  advantages  and  disadvantages. 
The  first  approach,  which  we  term  here  the  Erpenbeck 
equation  of  state  method  (E-EOS),  is  the  most  indirect 
of  the  approaches.  The  original  version  of  the  method 
involves  performing  several  separate  simulations  at 
appropriately  chosen  temperatures  and  pressures.  Each 
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simulation  generates  an  equation  of  state  point  for 
subsequent  evaluation  of  the  Hugoniot  expression 
followed  by  interpolation  to  locate  the  point  at  which 
the  expression  is  zero.  The  molecular  dynamics  (MD) 
method  has  been  implemented  in  the  Erpenbeck 
approach  using  reactive  potentials  that  mimic  chemical 
bond  breaking  and  forming  between  species  [12,  13]. 

An  approach  for  calculating  the  shock  Hugoniot 
properties  from  classical  molecular  simulation  that  is 
more  direct  than  the  E-EOS  method  is  the  piston-driven 
molecular  dynamics  method  [7,  8].  Piston-driven  MD 
generates  a  point  on  the  Hugoniot  curve  from  a  single 
simulation,  thus  avoiding  the  need  to  calculate  several 
EOS  points  (as  in  the  E-EOS  approach)  in  order  to 
obtain  the  desired  result.  The  method  mimics  the 
laboratory  system  by  calculating  properties  behind  the 
shock  discontinuity  in  a  shock  wave  simulation.  Shock 
waves  are  produced  by  hitting  the  free  edge  of  the 
molecular  solid  with  a  rigid  layer  of  atoms  that  are 
moving  at  a  constant  velocity.  Different  shocked  states 
are  obtained  by  starting  with  different  initial  piston 
velocities  [8]. 

A  third  approach,  termed  the  uniaxial  Hugoniostat 
method  [9],  is  a  molecular  dynamics  method  which 
utilizes  modified  equations  of  motion  that  constrain  the 
system  to  states  that  correspond  to  points  on  the  shock 
Hugoniot  curve.  This  method  is  computationally  more 
efficient  than  both  MD  implementations  of  the  E-EOS 
method  and  the  piston-driven  shock  wave  simulations. 
The  E-EOS  method  requires  a  system  of  only  a  few 
hundred  atoms  (with  periodic  boundary  conditions 
imposed),  but  several  simulations  are  required  to 
generate  a  single  shock  Hugoniot  point.  The  piston- 
driven  MD  method  will  produce  results  from  a  single 
simulation.  However,  the  system  size  must  be  suffi¬ 
ciently  large  so  that  the  properties  behind  the  propagat¬ 
ing  shock  wave  can  be  averaged  over  time.  The  uniaxial 
Hugoniostat  method,  on  the  other  hand,  can  generate  a 
point  on  the  Hugoniot  curve  from  a  single  simulation 
whose  system  size  is  relatively  small. 

Unfortunately,  none  of  these  MD  methods  can  be 
applied  to  the  calculation  of  the  shock  Hugoniot  locus 
over  a  wide  range  of  temperatures  and  pressures  unless 
certain  conditions  are  fulfilled.  Most  energetic  materials 
respond  to  shock  by  decomposing  into  a  complex 
(sometimes  heterogeneous)  mixture  of  many  different 
chemical  species.  Thus,  for  these  multi-component  sys¬ 
tems,  the  implementation  of  the  uniaxial  Hugoniostat 
method,  the  MD  E-EOS  method  or  the  piston-driven 
MD  method  requires  either:  (1)  a  priori  knowledge  of  the 
relative  concentrations  of  each  chemical  species  in  the 
shocked  state;  or  (2)  a  reactive  potential  that  simulates 
bond  breaking  and  bond  formation.  Typically,  the 
relative  species  concentrations  of  the  shocked  state  are 


lacking;  moreover,  knowledge  of  these  quantities  is 
desired.  Furthermore,  although  significant  advances 
have  been  made  in  developing  reactive  potentials  for 
shocked  materials,  the  potentials  are  presently  limited 
to  idealized  representations  of  the  chemistry  that  occurs 
[13-26].  The  most  accurate  interaction  potentials 
for  energetic  materials  available  at  this  time  are  all 
non-reactive  [27^11]. 

Monte  Carlo  methodologies  circumvent  some  of 
the  restrictions  associated  with  the  MD  methods  for 
calculating  shock  Hugoniot  states.  The  reactive  Monte 
Carlo  method  (RxMC)  [10]  and  the  composite  Monte 
Carlo  method  [11]  have  both  been  used  to  calculate 
Hugoniot  properties  through  the  E-EOS  approach. 
These  two  closely  related  methods  do  not  require  a 
reactive  potential  or  a  priori  knowledge  of  species 
concentrations  for  each  Hugoniot  state.  They  also  do 
not  require  the  specification  of  species  chemical  poten¬ 
tials  or  chemical  potential  differences  to  determine 
chemical  equilibrium  states  of  the  reactive  mixtures. 
Both  methods  have  been  applied  to  simple,  spherically- 
averaged  intermolecular  potentials  [10,  11]  but  can 
readily  be  applied  to  complex  potentials  that  include 
multi-site  and/or  electrically-charged  species  as  well 
as  multi-phase  mixtures.  Therefore,  in  the  absence 
of  reactive  potentials  or  a  priori  knowledge  of  species 
concentrations,  the  only  applicable  approaches  for 
simulating  the  Hugoniot  properties  of  a  shocked 
material  are  the  Erpenbeck  EOS  method  performed 
using  either  the  RxMC  or  composite  MC  methods. 

As  previously  mentioned,  however,  the  original  E-EOS 
method  requires  simulations  of  several  equation  of  state 
points  to  generate  a  single  point  on  the  shock  Hugoniot 
curve.  Each  separate  simulation  requires  sufficient 
equilibration  and  data  collection  steps.  In  an  effort  to 
reduce  the  number  of  steps  and  to  minimize  associated 
computational  costs,  we  have  implemented  a  numerical 
approach  within  the  framework  of  the  E-EOS  approach. 
The  resulting  method  requires  only  a  single  simulation 
to  determine  a  point  on  the  Hugoniot  curve.  The  fitting 
procedure  used  to  determine  the  root  of  the  Hugoniot 
expression  (Hg  =  0)  in  the  original  version  of  the 
E-EOS  approach  is  replaced  by  an  iterative  numerical 
procedure  built  into  the  framework  of  the  simulation. 

A  brief  illustration  of  the  method  using  isothermal- 
isobaric  ensemble  ( NPT )  simulations  is  given.  The 
simulation  is  initiated  at  the  specified  temperature  and 
pressure,  and  the  Hugoniot  expression  is  evaluated  (using 
instantaneous  values  of  V,  P  and  E  that  depend  only 
on  the  current  configuration)  at  periodic  intervals  during 
the  simulation  run  and  accumulated  for  averaging.  The 
averaged  Hugoniot  value  is  then  used  in  a  numerical 
root-finding  algorithm  (e.g.  Newton-Raphson  [42])  to 
provide  an  estimate  of  the  Hugoniot  pressure.  The 
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imposed  pressure  for  the  simulation  is  subsequently 
changed  to  correspond  to  the  new  estimate  of  the 
Hugoniot  pressure.  The  simulation  continues  using  this 
new  imposed  pressure  constraint.  This  process  is  repeated 
until  the  Hugoniot  function  converges  to  zero  (more 
precisely,  within  a  desired  tolerance  of  zero).  The  value  of 
the  pressure  and  corresponding  volume  averaged  over 
the  entire  simulation  run  characterize  the  Hugoniot  state 
at  that  temperature. 

The  method  is  akin  to  the  phase  equilibria  methods 
that  utilize  thermodynamic  integration  to  determine 
coexistence  behaviour  [43-45].  In  these  calculations,  a 
finite-difference  algorithm  is  used  to  numerically  inte¬ 
grate  the  differential  equations,  the  Clausius-Clapeyron 
[43,  44]  or  the  Gibbs-Duhem  expressions  [45],  which 
govern  the  changes  in  thermodynamic  parameters  along 
the  phase  coexistence  curve.  Similarly,  for  the  method 
introduced  here,  a  numerical  estimate  of  the  root  of  the 
Hugoniot  expression  is  made.  Both  approaches  are 
iterated  until  the  desired  convergence  has  been  reached. 

In  summary,  we  demonstrate  the  accuracy  and 
stability  of  a  method  to  calculate  the  shock  Hugoniot 
properties  of  materials  when  implementing  the  E-EOS 
method  either  in  an  MD  or  RxMC  framework. 
(Implementation  of  the  modified  E-EOS  method  using 
the  composite  MC  method  is  analogous  to  the  RxMC 
method  and  will  not  be  demonstrated  here.)  The 
method,  which  we  term  the  adaptive  Erpenbeck 
equation  of  state  method  (AE-EOS),  is  intended  to  be 
a  tool  that  is  an  alternative  to  the  existing  methodologies 
in  order  to  overcome  some  of  their  limitations.  We 
demonstrate  the  validity  of  the  AE-EOS  method  for 
calculating  the  Hugoniot  properties  of  liquid  N2  in  the 
non-reactive  regime  using  molecular  dynamics  and  in 
the  reactive  regime  using  reactive  Monte  Carlo.  The 
outline  of  the  paper  is  as  follows.  The  formalism  and 
practical  details  of  the  methodology  are  presented  in 
section  2.  Applications  of  the  method  are  given  in 
section  3,  while  assessments  of  the  results  and  possible 
extensions  of  the  method  are  given  in  section  4. 

2.  Methodology 

2.1.  Formalism 

The  Hugoniot  function  If  =  Hg[E,  T.  P( or  V)],  and 
so  a  search  for  If  =  0  in  the  original  Erpenbeck  EOS 
method  is  implemented  by  fixing  two  independent 
variables  in  a  simulation  and  calculating  the  remaining 
variable.  For  example,  when  simulating  with  an 
isothermal-isobaric  ensemble,  E  is  calculated  from  the 
simulation.  Typically,  several  simulations  must  be 
performed  using  various  choices  of  the  independent 
variables  to  generate  sufficient  points  so  that  the 
Hugoniot  state  can  be  obtained  through  interpolation. 
(Detailed  outlines  of  the  original  E-EOS  approach  using 


the  molecular  dynamics  and  the  RxMC  methods  can  be 
found  in  [12]  and  [10],  respectively.)  The  adaptive 
Erpenbeck  equation  of  state  method  presented  in  this 
work  eliminates  the  interpolation  procedure  by  using  a 
root-finding  algorithm  to  determine  Hg  —  0.  For  exam¬ 
ple,  the  expression  for  finding  the  root  of  a  function 
using  the  Newton-Raphson  method  is  [42] 


-Ch-  1 


f(Xn) 

/'(■*»)  ’ 


(2) 


where  f'(x„)  —  df(x„)/dx„,  f(xn)  =  Hg  and  we  choose 
either  xn  —  E,  T  or  P  (or  V).  The  basis  of  the  AE-EOS 
method  is  to  begin  a  simulation  at  an  initial  x„  and  after 
a  prescribed  number  of  simulation  steps,  the  quantities 
f(x„ )  and  f(x, ,)  are  calculated  and  xn+\  is  determined. 
The  simulation  continues  using  the  predicted  value  xn+\. 
This  procedure  is  repeated  until  the  results  converge. 

In  the  following,  we  demonstrate  the  AE-EOS 
method  for  P  as  the  independent  variable  used  to  find 
the  root  of  the  Hugoniot  function.  (Considerations 
for  choosing  the  other  variables  E,  T  and  V  as  the 
independent  quantity  in  the  root-finding  algorithm  are 
presented  in  the  Appendix.)  The  working  expression 
for  predicting  the  pressure  of  the  Hugoniot  state  within 
the  framework  of  the  Newton-Raphson  procedure 
(equation  (2))  is  then 


^predicted  _ ^current 


^current 

dHfneai/dP' 


(3) 


In  the  procedure  presented  here,  the  predicted  pressure, 
^predicted ^  js  determined  using  averaged  instantaneous 
values  of  the  Hugoniot  expression  (equation  (1))  and 
its  derivative  with  respect  to  pressure,  i.e.  HgUrrent  =  (Hg) 
and  d//gUrrent/d/)  =  {dHg/dP),  where  ‘O’  denotes 
averages  of  instantaneous  values.  The  instantaneous 
Hugoniot  values  (and  derivatives)  were  determined 
using  the  corresponding  instantaneous  values  of  P,  E 
and  T  generated  during  the  simulation. 

Next,  since  the  internal  energy  can  be  written  as  [10] 

9 

e=J2  y<H° + c/C0nf  -RT>  (4) 

i=  1 

then  equation  (1)  can  be  rewritten  as 


-  \(PVa  -  PV  +  P0  V0  -  VP0),  (5) 

where  E0,  P0  and  V0  are  the  values  of  the  initial  state 
(and  thus  constant)  and  E/conf  =  ]T  ]T  Uy(ry),  where  Uy 

i  j>i 


3312 


J.  K.  Brennan  and  B.  M.  Rice 


is  the  pair  potential  energy  [46].  All  quantities 
in  equations  (4)  and  (5)  are  used  here  in  the  context 
of  instantaneous  quantities  that  depend  only  on  the 
current  configuration.  The  derivative  term  required  in 
equation  (3)  is  then 


dH, 

d  P 


V 


i=i 


d  P' 


d  P 


1  d 

2cLP 


(PV0-PV  +  P0V0-VP0). 


(6) 


The  hrst  and  third  terms  on  the  right-hand  side  of 
equation  (6)  are  not  functions  of  P  and  can  be 
eliminated.  Further  since  C/cont  is  calculated  as  an 
instantaneous  value,  (7conf  =  Uconf(r)  only  and  thus  can 
be  eliminated.  Finally,  the  last  term  is  readily  solved,  so 
that  equation  (6)  reduces  to 


The  algorithm  for  the  adaptive  Erpenbeck  equation  of 
state  method  with  P  chosen  as  the  independent  variable 
is  as  follows. 


Step 

Step 

Step 


1:  Set  the  temperature  for  the  Hugoniot  state 
WHt). 

2:  Guess  the  pressure  for  this  Hugoniot  state 

^^Dcurrent^ 


3:  Perform  an  isothermal-isobaric  ensemble 
simulation  (MD  or  RxMC)  at  THg  and 

jDCurrent 


Step  4:  After  allowing  the  system  to  relax  to  Rcurrent, 
accumulate  instantaneous  values  of  Hg  and 
dHg/dP  during  the  simulation  using 
equations  (5)  and  (7),  respectively. 

Step  5:  After  a  prescribed  number  of  steps,  cal¬ 
culate  averaged  values  of  Hg  and  dHg/dP 
and  predict  the  Hugoniot  pressure  using 
equation  (3).  (The  averaged  values  of  Hg 
and  dHg/dP  can  be  determined  by  several 
methods,  which  are  considered  in  the  next 
section.) 

Step  6:  Repeat  steps  (3)-(6)  until  the  results  con¬ 
verge  to  the  desired  statistical  uncertainty. 


2.2.  Practical  details 

Next,  we  consider  a  few  practical  details  of  implement¬ 
ing  the  AE-EOS  method.  Our  intent  is  to  generalize  the 
method  for  implementation  into  any  of  the  standard 
molecular  simulation  techniques  (MD,  RxMC  or  com¬ 
posite  MC).  We  consider  the  effect  of  several  parameters 
on  the  accuracy  and  stability  of  the  method.  Below  we 
address  these  issues,  the  logic  behind  our  choices  and  the 
trade-offs  involved. 


2.2.1.  Root-finding  algorithm 

There  exists  a  wide  range  of  root-finding  algorithms, 
including  the  Newton-Raphson  method,  the  secant 
method,  the  bi-section  method  and  Halley’s  method 
[42].  Newton-Raphson  (equation  (2))  is  a  rather  straight¬ 
forward  method  but  can  be  unstable  near  a  horizontal 
asymptote  or  local  minimum.  A  similar  algorithm  is 
Halley’s  method  which  includes  an  additional  term 
(d/'(x„)/d.v„)  from  the  Taylor  series  in  the  derivation 
of  the  method.  When  the  pressure  is  chosen  as  the 
independent  variable,  d2Hg/dP 2  =  0  (see  equation  (7))  so 
Halley’s  method  reduces  to  the  Newton-Raphson 
method.  Another  root-finding  method  is  the  secant 
method,  which  estimates  the  derivative  term  using 
f'(xn)  =  [f(x„)-f(x„- 1  )]/(.\-„  -  x„-i).  However,  use  of 
the  method  in  the  AE-EOS  method  requires  two 
recent  points  along  the  Hugoniot  curve  as  opposed 
to  only  one  for  the  Newton-Raphson  method.  Finally, 
if  we  can  be  certain  that  the  solution  of  the  Hugoniot 
expression  lies  within  a  known  interval,  then  we  can 
iteratively  converge  to  the  solution  using  the  bi-section 
method.  However,  a  balance  must  be  established  between 
statistical  uncertainty  and  the  desired  convergence 
when  implementing  the  bi-section  method  in  a  molecular 
simulation,  since  statistical  fluctuations  cannot  be  greater 
than  the  size  of  the  interval.  In  an  effort  to  keep  the 
method  proposed  here  as  general  and  straightforward  as 
possible,  we  have  implemented  the  Newton-Raphson 
method.  The  well-known  problems  of  this  method  near  a 
local  minimum  or  asymptote  have  not  been  encountered 
for  the  Hugoniot  expression  in  this  work  as  well  as  for 
other  work  [10,  12,  13],  but  one  should  be  mindful  of  its 
limitations. 

2.2.2.  Initial  guess  of  pressure 

In  any  molecular  simulation,  it  is  necessary  to  design 
a  starting  configuration  so  that  the  relaxed  system  is 
physically  reasonable  and  computationally  consistent. 
For  example,  consider  an  isothermal-isobaric  ensemble 
simulation  where  periodic  boundaries  are  imposed 
and  where  the  potential  energy  function  has  a  limited 
interaction  range.  In  such  a  case,  an  appropriate  number 
of  molecules  must  be  chosen  so  that  the  relaxed  box  size 
is  consistent  with  the  minimum  image  convention,  i.e. 
one-half  the  box  size  must  be  greater  than  or  equal  to 
the  potential  cut-off  distance  [46].  Similarly,  appropriate 
starting  conditions  for  the  AE-EOS  method  are 
required,  particularly  for  the  initial  guess  of  the  imposed 
pressure.  Since  the  converged  result  (i.e.  Hg  —  0)  will 
produce  a  Hugoniot  pressure  that  is  equal  to  the 
imposed  pressure  (within  some  specified  tolerance),  it 
is  desirable  to  choose  an  initial  pressure  that  is  a  good 
estimate  of  the  actual  Hugoniot  pressure.  Although  we 
will  show  in  section  3  that  the  AE-EOS  method  is  largely 
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Table  1.  Estimates  of  the  initial  pressure  for  shocked  liquid  N2." 
Mole  fraction 


77  K 

x  (N2) 

x  (N ) 

F/cm3g  1 

E[ kJg"1 

^initial  guess/GP^ 

^Hg/OPa 

%  difference" 

536.2 

1.0 

0.0 

0.866 

0.09181 

2.82 

2.96 

4.7 

883.9 

1.0 

0.0 

0.866 

0.3762 

4.35 

4.74 

8.2 

2008.4 

1.0 

0.0 

0.866 

1.435 

10.1 

10.1 

0.0 

3912.4 

1.0 

0.0 

0.866 

3.385 

20.6 

18.1 

-13.8 

6778.1 

1.0 

0.0 

0.866 

6.434 

37.0 

29.9 

-23.7 

7963.0 

1.0 

0.0 

0.866 

7.745 

44.0 

36.0 

-22.2 

9557.7 

1.0 

0.0 

0.866 

9.622 

54.1 

47.0 

-15.1 

10185.4 

1.0 

0.0 

0.866 

10.41 

58.4 

52.6 

-11.0 

10935.2 

1.0 

0.0 

0.866 

11.4 

63.7 

60.4 

-5.5 

12588.9 

1.0 

0.0 

0.866 

13.8 

76.4 

81.1 

5.8 

"Quantities  shown  are  based  on  the  estimating  scheme  described  in  the  text. 
*Taken  from  [10]. 

f%  difference  =  100  x  (PHg-Amtiai  gu eSS)/P//g- 


insensitive  to  the  initial  pressure  guess,  extremely  poor 
initial  guesses  could  result  in  numerical  failures.  For 
example,  consider  the  simulation  of  a  Hugoniot  state 
such  that  T Hg  >  T0  (recall  that  T is  the  temperature  of 
the  chosen  Hugoniot  state  and  T0  is  the  temperature 
of  the  unshocked  material).  If  the  initial  guess  of  the 
pressure,  ^initial guess,  causes  an  expansion  of  the  simula¬ 
tion  cell  such  that  the  specific  volume  is  larger  than  the 
specific  volume  of  the  unshocked  material,  a  negative 
pressure  value  will  be  predicted.  From  a  practical 
standpoint,  this  is  an  unphysical  occurrence  since  it 
implies  that  the  material  has  expanded  upon  shock 
rather  than  being  compressed.  Furthermore,  from  a 
computational  standpoint,  the  simulation  cell  will  never 
converge  to  a  negative  imposed  pressure.  Such  an 
occurrence,  however,  is  analogous  to  choosing  a  starting 
configuration  that  relaxes  to  a  physically  unreasonable 
and  computationally  inconsistent  state. 

Consider  the  following  ad  hoc  approach  to  choosing  a 
reasonable  initial  guess  of  the  pressure  for  equation  (3). 
First,  assume  that  the  shocked  material  does  not 
decompose  (i.e.  chemically  react).  This  is  a  reasonable 
approximation  at  low  shock  pressures  and  reduces  the 
first  term  on  the  right-hand  side  of  equation  (5)  to 
^starting material-  NexF  neglect  the  contribution  of  C/conf 
and  thus  eliminate  the  second  term  on  the  right-hand 
side  of  equation  (5).  This  approximation  has  no  physical 
justification,  however  a  short  simulation  could  be 
performed  to  calculate  t/conf  (although  probably  unne¬ 
cessary  given  the  lack  of  sensitivity  of  the  final  result  on 
the  initial  guess  of  the  pressure).  Finally,  estimate  the 
amount  of  compression  the  starting  material  will 
undergo,  e.g.  V  —  0.7  VQ.  This  estimate  presumably  can 
be  predicated  on  previous  studies  of  the  material  or 
similar  materials.  With  these  approximations,  an  initial 
estimate  of  the  Hugoniot  pressure  can  be  determined. 


Estimates  of  the  initial  pressures  for  shocked  liquid 
N2  simulations  using  this  scheme  are  compared  with 
Hugoniot  pressures  calculated  in  an  earlier  study  [10]  in 
Table  1.  Small  differences  between  Pinitiai  guess  and  PHg 
are  found  at  low  pressures.  Obvious  improvements  at 
higher  pressures  can  be  made  by  increasing  the  estimate 
of  the  amount  of  compression  of  the  starting  material. 
More  importantly  however,  in  section  3,  we  assess  the 
effect  of  the  initial  guess  of  the  pressure  on  the  stability 
and  accuracy  of  the  AE-EOS  method.  As  will  be  shown, 
the  estimates  in  Table  1  provide  more  than  sufficient 
initial  guesses,  since  guesses  as  large  as  ±67%  Png 
converge  properly.  Furthermore,  as  points  along  the 
Hugoniot  curve  are  determined,  better  estimates  of  the 
initial  pressure  can  be  made  by  using  these  Hugoniot 
states. 

2.2.3.  Averaging 

For  completeness  of  study,  we  consider  three  different 
schemes  for  averaging  the  instantaneous  values  of  Hg 
and  dHg/dP ,  which  are  then  used  in  equation  (3):  (1) 
block  averages;  (2)  running  averages;  and  (3)  block- 
to-running  averages.  Block  averages  are  taken  from  a 
limited  number  of  configurations  immediately  preceding 
the  pressure  adjustment  step,  while  running  averages 
are  taken  continuously  over  all  configurations  generated 
during  the  simulation  run.  The  block-to-running 
averages  scheme  uses  a  block-averaging  scheme  for  the 
equilibration  period  of  the  simulation  run  and  then 
continues  with  a  running  average  scheme  for  the 
production  period.  This  scheme  may  most  effectively 
remove  the  effects  of  a  poor  initial  guess,  while  we 
expect  the  running  average  scheme  to  be  the  most 
effective  alternative  since  fluctuations  in  the  pressure 
will  become  increasingly  damped  as  the  simulation 
proceeds.  Block  averaging  methods  will  likely  be  more 


3314 


J.  K.  Brennan  and  B.  M.  Rice 


slowly  converging  at  best,  and  unstable  at  worst. 
Moreover,  running  average  schemes  have  been  the 
most  successful  scheme  in  the  finite-difference  algo¬ 
rithms  used  in  the  phase  coexistence  methods  mentioned 
previously  [43^45]. 

2.2.4.  Pressure  prediction  frequency 

We  also  consider  the  effect  of  the  frequency  of 
re-setting  the  pressure  during  the  simulation.  Less 
frequent  updates  are  expected  to  cause  the  results  to 
converge  more  slowly  while  more  frequent  updates 
could  possibly  cause  the  root-finding  scheme  to  become 
unstable  or  to  fluctuate  too  greatly. 

A  final  note  on  the  convergence  of  the  system  to 
the  predicted  pressure  is  considered.  Step  (4)  in  the 
algorithm  outline  (see  section  2.1)  allows  the  system  to 
converge  to  the  predicted  pressure  value  (within  a  few  % 
of  the  predicted  pressure  for  the  most  recent  simulation 
steps)  before  re-evaluating  the  Hugoniot  expression  and 
its  derivative  (dHg/dP)  in  the  equilibration  period  only. 
This  ensures  that  even  for  large  changes  in  the  predicted 
pressure,  equilibrated  information  is  still  used  in  the  //g 
and  dPfg/dP  calculation.  Typically,  these  large  changes 
will  only  occur  during  the  earliest  stages  of  the 
simulation.  At  later  times  during  the  production  cycles, 
this  criterion  is  nearly  always  satisfied. 


3.  Application:  shock  Hugoniot  states  of  liquid  N2 

For  demonstrative  purposes,  several  shock  Hugoniot 
states  of  liquid  nitrogen  are  considered.  The  shock 
Hugoniot  properties  were  predicted  based  on  the  initial 
states  calculated  previously:  T  —  77.0  K;  p  =  0.808 
gcm“3;  P=  50.49  MPa;  E  —  —0.441  kJ  g”1  [10],  At 
pressures  higher  than  ~30GPa  along  the  Hugoniot 
curve,  the  dissociation  reaction  of  molecular  nitrogen 
(N24^2N)  occurs.  Therefore,  we  demonstrate  the 
AE-EOS  method  using  the  molecular  dynamics  techni¬ 
que  only  at  pressures  below  30  GPa  while  we  demonstrate 
the  AE-EOS  method  using  the  RxMC  method  at  a 
wider  range  of  pressures.  Particles  interact  through  the 
exponential-six  potential,  which  can  be  expressed  as 


U exp— 6(0 


OO, 


1  —  (6/a) 


6 

-exp I  aj 


r  <  r core? 

r  >  r core' 


(8) 


where  e  is  the  depth  of  the  attractive  well  between 
particles,  rm  is  the  radial  distance  at  which  the  potential 
is  a  minimum,  while  a  controls  the  steepness  of  the 
repulsive  interaction.  The  cut-off  distance  rcore  is 
included  to  avoid  the  unphysical  singularity  in  the 
potential  function  as  r  ->  0.  rcore  is  the  smallest 


Table  2.  Exponential-6  potential  parameters  [47]. 


Species 

Y  core/  A 

^ra/A 

e/khk  1 

a 

n2 

1.13 

4.2005 

101.10 

12.684 

N 

0.98 

2.5688 

88.181 

11.013 

positive  value  for  which  dUexp-6(r)/dr  =  0  and  is 
obtained  by  iterative  solution  of  equation  (8).  The 
potential  parameters  for  the  species  considered  in  this 
work  are  given  in  Table  2.  A  spherical  cut-off  for 
the  particle-particle  interactions  was  applied  at  2.5rm  N2 
with  long-range  corrections  added  to  account  for 
interactions  beyond  this  distance  [48].  Electrostatic 
interactions  between  species  were  ignored.  The  unlike 
interactions  between  species  i  and  j  were  approximated 
by  the  Lorentz-Berthelot  combining  rules  [49]  for  Ey,  ay 
and  rnvy: 

£//  —  (£,£/)  ;  ay  —  (a,aj)  ;  rm  y  =  (r my  -}-  rmy) /'l ; 

(9) 


while 

A  ore,//  —  ( /'core./  T  C  core\/)/2.  (10) 

3375  N?  molecules  were  used  with  all  calculated  quan¬ 
tities  reduced  by  the  exponential-6  potential  energy  (e) 
and  size  (rm)  parameters  of  N2.  Periodic  boundary 
conditions  were  imposed  for  all  dimensions.  Thermo¬ 
chemical  reference  data  were  used  in  calculating 
the  ideal-gas  enthalpies  (H°)  required  in  equation  (5) 
[50,  51], 


3.1.  Molecular  dynamics 
3.1.1.  Simulation  details 

Molecular  dynamics  simulations  in  the  isothermal- 
isobaric  ensemble  were  performed  using  the  leap-frog 
Verlet  algorithm  [46,  48]  and  the  Melchionna  modifica¬ 
tion  of  the  Hoover-Nose  equations  of  motion  [52],  A 
thermostatting  rate  of  50  ps-1  was  used  to  maintain  the 
imposed  temperature  while  a  barostatting  rate  ranging 
from  0.032-0.042  ps-1  was  used  to  maintain  the  imposed 
pressure,  initial  configurations  were  generated  from  a 
face-centred-cubic  (fee)  lattice  structure  with  initial 
particle  velocities  selected  from  a  Boltzmann  distribu¬ 
tion  that  corresponded  to  the  imposed  temperature. 
Preceded  by  an  equilibration  period  of  0.127-0.254  ns 
during  which  the  pressure  was  not  re-set  using 
equation  (3),  trajectories  were  followed  for  1.32  ns 
with  time  steps  ranging  from  0.00763-0.0102  ps.  All 
pressure  values  reported  were  determined  using  the  virial 
theorem  [46], 
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3.1.2.  AE-EOS  method  details 

Three  state  points  along  the  shock  Hugoniot  curve 
were  determined:  T  =  883.9;  3912.4;  6778.1  K.  These 
state  points  are  below  the  regime  in  which  N2  dissociates 
into  atomic  nitrogen.  For  each  point,  the  Hugoniot 
pressure  (P11A  was  predicted  in  two  simulations,  one  in 
which  the  initial  pressure  was  much  lower  than  the 
Hugoniot  pressure  and  one  in  which  the  initial  pressure 
was  too  high.  The  values  of  the  initial  guesses  are  given 
in  Table  3  and  correspond  to  differences  in  the  known 
value  [10]  by  ±67%.  The  effect  of  the  frequency  of 
re-setting  the  imposed  pressure  was  also  studied.  Two  cases 
were  considered,  re-setting  at:  (a)  every  100  steps;  and 
(b)  every  500  steps.  Following  the  initial  equilibration 
period  used  to  relax  the  system  from  the  fee  crystal  to 
the  imposed  thermodynamic  condition,  an  additional 
0.305  ns  of  the  trajectory  was  used  to  further  equilibrate 
the  system  after  the  AE-EOS  algorithm  is  implemented 
(i.e.  the  pressure  is  re-set  at  specified  intervals).  All 
quantities  calculated  during  this  time  interval  were  not 
included  in  the  final  averages.  A  tolerance  value  of  ±5% 
was  used  in  Step  (4)  for  the  pressure  (see  section  2.1),  i.e. 
the  calculated  pressure  was  required  to  be  within  ±5% 
of  the  most  recent  PHjf  prediction  before  re-evaluating 
Hg  and  dHg/dP  and  re-setting  of  the  imposed  pressure. 

3.1.3.  Results 

A  comparison  between  the  Hugoniot  properties 
predicted  using  the  original  E-EOS  and  AE-EOS 
methods  is  shown  in  table  3.  Good  agreement  is  found 
for  all  cases  considered,  with  pressure  and  specific 


volume  values  well  within  statistical  uncertainty.  Figures 
1(a)  and  (b)  are  the  instantaneous  values  every  100th 
time-step  for  the  Hugoniot  expression  (equation  (5)) 
and  the  predicted  pressure  (equation  (3)),  respectively,  at 
the  Hugoniot  point  T  —  883.9  K.  The  results  converge 
within  a  few  thousand  time-steps.  The  figures  show  that 
Hg  fluctuates  about  zero  and  the  predicted  pressures 
converge  to  PHg. 

Table  3  also  shows  that  no  dependence  on  the  initial 
pressure  guess  is  found,  with  values  within  a  few  percent 
of  the  E-EOS  method  results.  A  graphical  demonstra¬ 
tion  of  the  stability  of  the  algorithm  is  given  in 
Figures  1  (c)  and  (d)  for  the  T  =  6778.1  K  point  along 
the  Hg  curve.  In  figure  1  (c)  the  fluctuations  of  the  Hg 
expression  are  shown  to  increase  sporadically  through¬ 
out  the  simulation  but  the  scheme  ‘recovers’  within  a  few 
thousand  time-steps.  A  plot  of  the  predicted  pressure  for 
this  point  (figure  1  (d))  shows  that  the  pressure  predicted 
in  the  Newton-Raphson  scheme  remains  stable  during 
these  occasional  fluctuations  of  the  Hg  value. 

As  expected,  fluctuations  are  more  balanced  when 
the  Hg  expression  is  evaluated  less  frequently  since 
the  system  configuration  is  allowed  longer  to  converge 
to  the  predicted  pressure.  This  is  evident  by  comparing 
figures  1  (c)  and  (e),  where  the  difference  between 
the  data  shown  is  based  on  the  frequency  of  the  Hg 
evaluation  and  imposed  pressure  re-set:  figure  1  (c)  is 
for  every  100  time-steps,  while  figure  1  (e)  is  for  every 
500  time-steps.  While  fluctuations  may  be  steadier  for 
less  frequent  evaluations,  still  no  noticeable  dependen¬ 
cies  of  the  final  values  of  Pnt  and  Vng  on  the  frequency 


Table  3.  Predicted  shock  Hugoniot  states  of  liquid  N2  using  molecular  dynamics  in  the  AE-EOS  method." 

Hg  evaluation  frequency 


Every  100  steps  Every  500  steps 


Anitiai/  F/cm3/  V  /cm3 

GPa  r/ K  P/GPa  mol  N2  Hf  kjg’1  Tj K  P/GPa  mol-1  N2  #g/kJg_1 


hT  =  883.9  K  ;  P  =  4.74  GPa  V=  19.82  cm3  mole-1  N2 


1.56 

883.9  (1.2) 

4.82  (9) 

19.7  (2)  — 2.356E-4  (0.013) 

883.9  (1.1) 

4.8  (1) 

19.9  (3) 

8.901E-3 

(0.16) 

7.92 

883.9  (1.2) 

4.81  (2) 

19.74  (3)  1 .001 E-4  (0.023) 

883.9  (1.2) 

4.8(1) 

20.0  (3) 

1.600E-2 

(0.19) 

bT=  3912.4  K  ;  P=  18.1  GPa  V= 

15.57  cm3  mole-1  N2 

5.97 

3912  (52) 

18.4  (3) 

15.6  (2)  -1.550E-3  (0.15) 

3912  (52) 

18.4  (2) 

15.6  (2) 

— 3.101E-5 

(0.043) 

30.23 

3912  (53) 

18.4  (1) 

15.56  (8)  2.790E-3  (0.063) 

3912  (54) 

17.8  (6) 

16.0  (4) 

1.841E-1 

(1.1) 

hT=  6778.1  K  ;  P  =  29.9  GPa  V= 

14.05  cm3  mole-1  N2 

9.87 

6778  (94) 

29.9  (4) 

14.1  (2)  — 1.867E-3  (0.25) 

6778  (94) 

29.9  (9) 

14.1  (4) 

1.511E-1 

(1.6) 

49.93 

6778  (93) 

29.9  (5) 

14.1  (2)  — 5.581E-3  (0.35) 

6778  (94) 

29.9  (5) 

14.1  (2) 

1.355E-4 

(0.13) 

"Quantities  are  ensemble  averages.  Uncertainties  in  units  of  the  last  decimal  digit  are  given  in  parentheses:  e.g.,  883.9  (1.2)  means 
883.9  ±  1.2,  except  for  the  uncertainties  of  Hg  where  the  values  given  in  parentheses  are  absolute  values.  Uncertainties  reported  were 
determined  front  the  standard  deviation  of  the  instantaneous  values. 

^Hugoniot  states  taken  from  [10]. 
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Figure  1.  Instantaneous  simulation  results  for  predicted  points  along  the  Hugoniot  curve:  ( a)-(b )  T  =  883.9  K  point  with  a 
frequency  of  Hg  evaluation  and  subsequent  pressure  re-set  every  100  time  steps  and  initial  pressure  guess  of  1.56  GPa;  (c)-(d) 
T  =  6778.1  K  point  with  a  frequency  of  Hg  evaluation  and  subsequent  pressure  re-set  every  100  time  steps  and  initial  pressure 
guess  of  49.93  GPa;  same  as  (c)-(rf),  except  the  frequency  of  Hg  evaluation  and  subsequent  pressure  re-set  is  every  500 

time  steps.  The  left-most  frame  of  each  set  shows  the  value  of  the  Hugoniot  expression  calculated  from  equation  (5)  and  the 
right-most  frame  is  the  pressure  predicted  by  equation  (3).  Inset  in  (e)  shows  fluctuations  of  Hg  while  the  inset  in  (d)  shows  the 
fluctuations  in  the  pressure  for  an  isothermal-isobaric  ensemble  molecular  dynamics  simulation  carried  out  at  P  =  29.9  GPa 
and  T=  6778.1  K. 
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of  the  Hg  evaluation  and  subsequent  re-setting  of  the 
pressure  are  found.  For  comparison,  the  inset  in 
figure  1  ( d )  shows  the  fluctuations  in  the  calculated 
pressure  for  an  isothermal-isobaric  ensemble  molecular 
dynamics  simulation  carried  out  at  P—  29.9  GPa  and 
T  =  6778.1  K. 


3.2.  Reactive  Monte  Carlo 

3.2.1.  Simulation  details 

The  reactive  Monte  Carlo  method  was  used  to 
assess  the  accuracy  of  the  AE-EOS  method  at  a  wider 
range  of  conditions  than  considered  using  the  molecular 
dynamics  technique  including  conditions  under  which 
N2  dissociates  (N2<£>2N).  Details  of  the  methodology 
can  be  found  in  the  original  papers  [53-55]  as  well  as  in 
recent  applications  of  the  technique,  which  implemented 
the  E-EOS  method  [10].  In  addition  to  intermolecular 
potentials  that  describe  non-reactive  interactions 
between  species  N2  and  N  in  the  equilibrium  mixture, 
RxMC  also  requires  inputting  the  ideal-gas  internal 
modes  (vibration,  rotation,  electronic).  The  vibrational 
and  rotational  contributions  to  the  ideal-gas  partition 
functions  of  both  species  were  calculated  using  a 
standard  source  [56],  and  supplemented  with  electronic 
level  constants  that  included  the  ground  state  and 
six  excited  electronic  states  for  N?  [50],  along  with 
electronic  energy  levels  for  N  taken  from  Moore  and 
Gallagher  [51]. 

Constant-pressure  RxMC  simulations  of  shocked  N2 
were  initiated  from  3375  N2  particles  placed  on  an  fee 
lattice  structure.  Simulations  were  performed  in  steps, 
where  a  step  (chosen  with  equal  probability)  was  either 
a  particle  displacement,  a  forward  reaction  step,  or  a 
reverse  reaction  step.  A  change  in  the  simulation  cell 
volume  was  attempted  every  500  steps.  Simulations 
were  equilibrated  for  1.5  xlO6  steps  after  which  aver¬ 
ages  of  the  quantities  were  taken  over  8  x  106  steps. 
Uncertainties  were  estimated  using  the  method  of  block 
averages  by  dividing  the  production  run  into  10  equal 
blocks  [48],  Reported  uncertainties  are  one  standard 
deviation  of  the  block  averages.  (Note  that  these  block 
averages  do  not  correspond  to  the  averaging  used  in  the 
Pu  prediction  scheme.  )  The  maximum  displacement  and 
volume  change  were  adjusted  to  achieve  an  acceptance 
fraction  of  approximately  0.33  and  0.5,  respectively. 
Depending  on  the  system  conditions,  the  acceptance 
fraction  of  the  reaction  steps  ranged  from  0.075-0.375. 

3.2.2.  AE-EOS  method  details 

Three  points  along  the  shock  Hugoniot  curve  were 
considered:  T  =  2008.4;  7963.03;  10935.2  K.  Again  for 
each  point,  Ph  was  predicted  from  two  initial  guesses 


(values  are  given  in  Tables  4-6)  that  differ  from  the 
known  value  [10]  by  ±67%.  We  predicted  the  shock 
Hugoniot  properties  based  on  the  same  initial  state  used 
previously  [10]  and  in  the  molecular  dynamics  study 
above.  A  tolerance  value  of  ±2.5%  was  used  in  Step  4 
for  the  pressure  (see  section  2.1),  a  value  slightly  lower 
than  used  in  the  molecular  dynamics  simulations 
above.  Analogous  to  the  study  above,  the  effect  of  the 
frequency  of  re-evaluating  the  Hugoniot  pressure  and 
subsequent  re-setting  of  the  imposed  pressure  was 
studied.  Two  cases  were  considered,  re-setting  at:  (a) 
every  5000  steps;  and  (b)  every  50000  steps.  Addition¬ 
ally,  three  different  averaging  schemes  of  the  instanta¬ 
neous  values  of  Pfe  and  dPlg/dP  were  assessed:  (a) 
block  averages;  (b)  running  averages;  and  (c)  block-to- 
running  averages.  A  description  of  each  type  is  given 
in  section  2.2.  For  each  of  the  initial  pressure  guesses, 
all  six  series  (2  Hg  frequency  evaluations  x  3  averaging 
schemes)  were  considered. 

3.2.3.  Results 

Comparisons  between  the  Hugoniot  properties  pre¬ 
dicted  by  the  E-EOS  and  the  AE-EOS  methods  are  given 
in  Tables  4-6.  We  find  excellent  agreement  between  the 
E-EOS  and  the  AE-EOS  results  for  all  quantities 
calculated.  Nearly  all  AE-EOS  quantities  fall  within 
±0.5%  of  the  E-EOS  calculations  with  none  greater 
than  ±0.9%.  Consequently,  no  dependence  on  the  initial 
pressure  guess  or  the  Hugoniot  expression  evaluation 
frequency  is  evident.  Likewise,  no  apparent  dependence 
on  each  of  the  averaging  schemes  is  found.  A  detailed 
examination  of  the  instantaneous  values  reveals  the 
behaviour  of  the  averaging  schemes.  Figure  2  shows 
the  instantaneous  values  of  the  Hugoniot  expression 
throughout  the  simulation  run  for  the  T  =  2008. 4  K 
point  using  an  initial  pressure  guess  of  16.9  GPa.  Note 
that  to  discern  the  data  more  clearly,  the  ‘running’  and 
‘block-to-running’  curves  are  offset  on  the  y  axis  by  0.05 
and  0.1  kJg-1,  respectively.  Updates  of  the  imposed 
pressure  every  50  000  steps  (figure  2  (a))  and  every  5000 
steps  (figure  2(b))  are  shown.  Figures  2(a)  and  (b) 
exhibit  similar  behaviour  with  the  expected  increased 
statistical  ‘noise’  for  figure  2(b).  The  block-averaging 
scheme  fluctuates  throughout  the  simulation  runs  but 
in  a  consistent  and  stable  fashion  while  the  running- 
average  scheme  fluctuates  minimally.  The  block-to- 
running  averaging  scheme  exhibits  a  combination  of 
these  behaviours,  balanced  and  steady  fluctuation 
during  the  block-averaging  portion  followed  by  limited 
fluctuation  during  the  running  portion.  Insets  in 
figures  2(a)  and  iff)  show  fluctuations  for  the  running 
average  scheme  on  a  reduced  scale.  Similar  fluctuations 
are  found  for  the  ‘running’  portion  of  the  block-to- 
running  average  scheme.  Other  series  points  exhibit 
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Table  4.  Predicted  shock  Hugoniot  states  of  liquid  N2  at  T=  2008.4  K  using  the  reactive  Monte  Carlo  method.  Values  determined 
previously  [10]  by  the  E-EOS  method  are  P—  10.1  GPa  ;  V=  17.41  cm3 mole-1  N2;  x  (N2)=  1.00.“ 

Hg  evaluation  frequency 

Every  5000  steps  Every  50000  steps 


initial/  V/  V/ 


GPa 

P/GPa 

cm3  mol  1  N2 

bx  (N2) 

HgfkJ  g-1 

P/  GPa 

cm3  mol  'N2 

-v  (N2) 

tfg/kjg-1 

Block 

3.3 

averages 

10.1  (1) 

17.35  (3) 

1.00  (0) 

2.177E-5  (0.013) 

10.1  (1) 

17.35  (6) 

1.00  (0) 

— 4.706E-6  (0.16) 

16.9 

10.1  (1) 

17.35  (5) 

1.00  (0) 

— 9.019E-7  (0.023) 

10.1  (1) 

17.35  (4) 

1.00  (0) 

— 2.479E-5  (0.19) 

Running  averages 
3.3  10.1  (1) 

17.32  (1) 

1.00  (0) 

— 2.716E-5  (0.15) 

10.2  (1) 

17.33  (TO) 

1.00  (0) 

—  1.357E-4  (0.043) 

16.9 

10.2  (1) 

17.32  (13) 

1.00  (0) 

— 2.549E-5  (0.063) 

10.2  (1) 

17.33  (9) 

1.00  (0) 

—  1.648E-4  (1.1) 

Block 

3.3 

— >■  Running  averages 

10.1  (1)"  17.35  (9) 

1.00  (0) 

8.119E-6  (0.25) 

10.1  (1) 

17.34  (11) 

1.00  (0) 

1.816E-5  (1.6) 

16.9 

10.1  (1) 

17.35  (8) 

1.00  (0) 

— 4.286E-6  (0.35) 

10.1  (1) 

17.35  (7) 

1.00  (0) 

—  1.508E-5  (0.13) 

"Quantities  are  ensemble  averages.  Reported  uncertainties  shown  in  parentheses  are  one  standard  deviation  of  the  block  averages 
[48]  and  are  given  in  units  of  the  last  decimal  digit,  e.g.,  17.35  (3)  means  17.35  ±0.03. 

*Mole  fraction  of  N2,  so  .v(N2)  =  VN2/iVtotai  and  x(N)=  l/2iVN/7Vtotai,  where  7Vtotai  =  3375. 


Table  5.  Predicted  shock  Hugoniot  states  of  liquid  N2  at  T=  7963.0  K  using  the  reactive  Monte  Carlo  method.  Values  determined 
previously  [10]  by  the  E-EOS  method  are  P  =  36.0GPa  ;  V=  13.35 cm3 mole-1  N2;  .v(N2)  =  0.975." 


Hg  evaluation  frequency 


-^initial/ 

GPa 

Every  5000  steps 

Every  50000  steps 

P/ GPa 

VI 

cm3  mol  1  N2 

-v  (N2) 

HgfkJ  g-1 

P/GPa 

VI 

cm3  mol  1  N2 

.v  (N2) 

HgfkJ  g-1 

Block  averages 

11.9 

36.1  (1) 

13.33  (6) 

0.975  (1) 

2.129E-5  (0.013) 

36.2  (1) 

13.32  (3) 

0.975  (1) 

— 5.957E-4 

(0.16) 

60.1 

36.1  (1) 

13.33  (6) 

0.975  (1) 

1.118E-4  (0.023) 

36.1  (1) 

13.33  (4) 

0.975  (1) 

2.937E-5 

(0.19) 

Running  averages 

11.9 

36.2  (1) 

13.32  (8) 

0.975  (1) 

— 2.435E-5  (0.15) 

36.1  (1) 

13.33  (6) 

0.975  (1) 

5.444E-6 

(0.043) 

60.1 

36.2  (1) 

13.32  (6) 

0.975  (1) 

—  1.022E-5  (0.063) 

36.1  (1) 

13.33  (6) 

0.975  (1) 

— 2.631E-5 

(1.1) 

Block  - 

*  Running  averages 

11.9 

36.1  (1) 

13.33  (4) 

0.975  (1) 

1.942E-5  (0.25) 

36.1  (1) 

13.33  (4) 

0.975  (1) 

— 6.005E-4 

(1.6) 

60.1 

36.1  (1) 

13.33  (2) 

0.975  (1) 

2.088E-5  (0.35) 

36.1  (1) 

13.33  (4) 

0.975  (1) 

—  1.805E-4 

(0.13) 

"See  table  4  for  details. 


similar  behaviour,  though  again  with  no  particular 
averaging  scheme  demonstrating  noticeable  results.  This 
behaviour  gives  testament  to  the  AE-EOS  method’s 
overall  stability  and  robustness. 


4.  Discussion 

We  have  presented  a  computationally  efficient  meth¬ 
odology  for  calculating  the  shock  Hugoniot  properties 
of  materials  using  classical  molecular  simulations.  The 
method  is  an  extension  of  the  Erpenbeck  EOS  approach 
and  allows  for  the  determination  of  a  point  along  the 


Hugoniot  curve  from  a  single  simulation.  The  method, 
termed  the  adaptive  Erpenbeck  equation  of  state 
method  (AE-EOS),  uses  a  numerical  estimate  of  the 
root  of  the  Hugoniot  expression  to  determine  the 
corresponding  thermodynamic  state.  AE-EOS  is  applic¬ 
able  for  any  simulation  method  that  determines  points 
along  the  Hugoniot  curve  by  generating  EOS  points, 
which  includes  molecular  dynamics  [46,  48],  reactive 
Monte  Carlo  [53,  54]  and  the  composite  Monte  Carlo 
method  [11],  The  method  was  demonstrated  to  be 
accurate  and  numerically  stable.  For  the  systems  in 
this  study,  the  method  was  not  particularly  sensitive 
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Table  6.  Predicted  shock  Hugoniot  states  of  liquid  N2  at  T=  10935.2  K  using  the  reactive  Monte  Carlo  method.  Values 
determined  previously  [10]  by  the  E-EOS  method  are  _P  =  60.4  GPa  ;  V=  10.97  cm3  mole-1  N2;  x(N2)  =  0.821." 

Hg  evaluation  frequency 

Every  5000  steps  Every  50000  steps 


initial/  VI  V/ 


GPa 

PI  GPa 

cm3/mol  1  N2 

v  (N2) 

Hg/kJg-1 

P/GPa 

cm3  mol  1  N2 

x  (N2) 

Hg/kJg-1 

Block 

20.0 

averages 

60^6  (1) 

10.95  (9) 

0.820  (2) 

9.934E-5  (0.013) 

60.7  (1) 

10.95  (9) 

0.820  (2) 

2.759E-3  (0.16) 

100.9 

60.6  (1) 

10.95  (1) 

0.820  (3) 

2.528E-5  (0.023) 

60.6  (1) 

10.96  (8) 

0.820  (1) 

— 8.182E-4  (0.19) 

Running  averages 
20.0  60.7  (1) 

10.95  (6) 

0.820  (2) 

— 5.219E-5  (0.15) 

60.6  (1) 

10.95  (4) 

0.819  (1) 

— 4.466E-5  (0.043) 

100.9 

60.7  (1) 

10.95  (3) 

0.820  (1) 

— 3.152E-5  (0.063) 

60.7  (1) 

10.94  (5) 

0.820  (2) 

— 4.217E-4  (1.1) 

Block 

20.0 

— >■  Running  averages 
60.5(1)"  10^96  (6) 

0.820  (1) 

2.017E-5  (0.25) 

60.6  (1) 

10.96  (7) 

0.820  (2) 

1.909E-3  (1.6) 

100.9 

60.7  (1) 

10.95  (7) 

0.820  (2) 

—  1.463E-4  (0.35) 

60.6  (1) 

10.95  (3) 

0.820  (1) 

2.541E-4  (0.13) 

"See  table  4  for  details. 


to  the  algorithm  parameters,  indicating  the  method’s 
robustness  and  ability  to  be  readily  implemented. 
Furthermore,  substantial  savings  in  the  computational 
requirements  are  gained  through  the  implementation 
of  the  AE-EOS  method.  For  comparable  statistical 
uncertainty,  four  to  five  complete  simulations  are 
required  to  generate  a  single  point  on  the  Hugoniot 
curve  for  the  original  Erpenbeck  method,  while  only  a 
single  simulation  is  needed  in  the  AE-EOS  method. 
Hence,  a  computational  gain  of  approximately  four- 
five-fold  is  made  when  implementing  the  AE-EOS 
method.  Overall,  the  AE-EOS  method  appears  to  be 
slightly  more  accurate  and  stable  when  applied  to  the 
RxMC  method  as  compared  to  the  molecular  dynamics 
method.  This  may  be  attributed  to  the  relative  stability 
of  the  Monte  Carlo  method  over  the  molecular 
dynamics  technique  for  constant-temperature  and  con¬ 
stant-pressure  simulations. 

The  AE-EOS  method  is  intended  as  an  alternative 
tool  to  similar  techniques,  namely  piston-driven  molec¬ 
ular  dynamics  [7,  8]  and  uniaxial  Hugoniostat  molecular 
dynamics  [9],  since  AE-EOS  offers  some  notable 
advantages  when  simulating  reactive  mixtures.  The 
most  striking  advantage  is  that  the  AE-EOS  method 
does  not  require  either:  (1)  a  priori  knowledge  of  the 
relative  concentrations  of  each  chemical  species  in  the 
shocked  state;  or  (2)  a  reactive  potential  that  simulates 
bond  breaking  and  bond  formation.  At  least  one  of 
these  conditions  must  be  met  for  the  other  methods  to 
be  applicable  to  the  simulation  of  reactive  mixtures. 

Extensions  as  well  as  different  adaptations  of  the  AE- 
EOS  method  are  possible.  For  example,  as  outlined  in 
section  2.1,  there  exist  several  choices  of  independent 
variables  that  the  ensemble  can  be  based  upon  (these 


choices  are  considered  in  the  Appendix).  Furthermore, 
as  discussed  in  section  2.2  other  numerical  root-finding 
algorithms  that  converge  the  simulation  to  Hg  =  0  are 
possible.  Simplicity  and  generality  in  implementing  the 
method,  however,  has  led  to  the  choices  presented  here. 
Besides  the  numerical-based  approaches  described  in 
this  work,  other  approaches  based  on  the  modification 
or  invention  of  a  statistical  ensemble  may  be  possible. 
For  example,  through  the  development  of  an  NVE- 
RxMC  ensemble,  where  E  —  E0  +  ^(P  +  P0)(V0  —  V); 
N  and  V  are  held  constant,  while  E  will  fluctuate  (the 
kinetic  energy  will  be  kept  nearly  constant).  In  addition 
to  the  standard  RxMC  moves,  attempts  to  randomly 
change  E  will  be  made  with  the  acceptance  of  the  move 
based  on 


P  cx 


for  Hg  <  A , 


(ID 


where  e  is  an  arbitrary  parameter  defining  the  width  of 
the  A  function.  The  method  would  be  similar  to  N VE  - 
MC  [46]  but  with  the  additional  RxMC-type  reaction 
moves.  A  still  more  rigorous  and  fundamental  approach 
is  the  development  of  a  constant-/7g  ensemble,  i.e. 
the  development  of  an  7VF//g-RxMC,  NPHg-RxMC,  or 
NTHg-RxMC  ensemble.  RxMC  moves  would  be 
accepted  based  on  criteria  derived  specific  to  these 
ensembles.  Although  in  principle  it  may  be  possible  to 
determine  an  expression  for  the  partition  functions 
for  such  ensembles,  it  may  not  be  possible  to  derive 
expressions  with  a  practical  form  (e.g.  which  neatly 
separate  inter-  and  intramolecular  terms).  If  such  a 
method  is  plausible,  it  would  be  analogous  to  the 
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Figure  2.  Instantaneous  values  of  the  Hugoniot  expression  for  the  T  =  2008.4  K  point  using  the  reactive  Monte  Carlo  method 
where  the  initial  pressure  guess  is  16.9  GPa.  Frequency  of  Hg  evaluation  is  shown  for:  (a)  every  50  000  steps;  and  (b)  every  5000 
steps.  Curves  corresponding  to  ‘running’  and  ‘block-to-running’  notation  in  the  legend  have  been  shifted  by  0.05  and  0.1  kJg~', 
respectively. 


uniaxial  Hugoniostat-MD  method.  It  remains  to  be 
seen,  however,  if  such  methods  would  offer  any 
computational  advantage  over  the  current  methods. 
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Appendix 

In  section  2.1  we  derived  the  AE-EOS  method  for 
pressure  as  the  independent  variable  in  the  root-finding 
algorithm.  Let  us  consider  the  other  possible  choices 
of  the  independent  variable  when  implementing  this 
approach,  namely,  E,  T  and  V.  The  choice  of  any  of 
these  variables  requires  an  ensemble  for  which  this 
particular  variable  is  fixed.  Therefore,  when  x„  —  E,  a 
constant-E  ensemble  is  required.  The  original  molecular 
dynamics  technique  is  based  upon  such  an  ensemble 
(NVE)  [46],  Furthermore,  a  direct  Monte  Carlo  method 
for  atomic  models  in  an  NVE  ensemble  has  been 
developed  [57,  58]  and  was  recently  extended  to 
molecular  models  that  include  all  internal  modes  in  the 
internal  energy  [59].  However,  the  choice  of  constrain¬ 
ing  E  is  somewhat  less  convenient  when  calculating 
points  along  the  Hugoniot,  since  typically  experimental 
measurements  provide  P-V/V0  or  P  T  curves. 

The  choice  of  x„  —  V  can  be  eliminated  due  to 
the  possibility  of  encountering  spurious  circumstances 
during  a  molecular  simulation.  For  example,  for  each 
iteration  of  the  root-finding  algorithm  a  new  simulation 
cell  volume  will  be  predicted.  However,  a  change  in  the 
simulation  cell  volume  during  a  molecular  simulation 
typically  requires  a  rescaling  of  the  particle  positions 
(techniques  exist  which  avoid  total  particle  rescaling  but 
apply  to  less  general  cases).  Therefore,  in  cases  where 
the  predicted  volume  is  considerably  smaller  than  the 
current  volume  (which  may  occur  during  the  first  step 
when  a  poor  initial  guess  of  the  Hugoniot  volume  is 
made)  it  is  likely  that  a  significant  number  of  particles 
will  have  substantial  overlap  and  create  an  energetically 
unfavourable  configuration.  For  Monte  Carlo  simula¬ 
tions,  this  is  less  prohibitive  but  for  molecular  dynamics 
simulations  such  a  situation  could  cause  significant 
deviations  in  the  conservation  of  energy  and  momen¬ 
tum  constraints.  For  this  reason  V  as  the  independent 
variable  in  the  AE-EOS  method  may  not  be  appropriate. 

Finally  for  x„  —  T,  consider  the  derivative  term 
required  in  equation  (2) 


since  77;°  can  be  written  as  a  function  of  T  (and  often  is 
determined  from  the  heat  capacity,  Cp,  in  this  form). 
The  second  term  on  the  right-hand  side  is  eliminated 
from  equation  (Al)  since  t/conf  =  Uconf(r)  only,  while 
the  third  term  reduces  to  R.  The  last  term  can  be 
expanded  so  that 


1  d 
2dE 


(PV0-PV  +  P0V0-  VPo) 


(A2) 


Simplifying  equation  (A2),  we  see  that  the  first  two 
terms  on  the  right-hand  side  require  further  evaluation 
of  a  dP/dT  term,  while  the  third  term  is  eliminated  since 
Po  and  Vo  are  constant.  If  we  consider  the  standard 
virial  equation  typically  used  in  molecular  simulation  to 
calculate  an  instantaneous  pressure  (see  e.g.  [46]),  then 


P  —  P  ideal  +  P e: 

gas 


,T  1^^  d  uconi 

—  pkT  —  -  2^  2_^  r~ 


i  j>i 


dr 


(A3) 


where  the  summations  are  taken  over  all  i-j  particle 
pairs.  Calculating  the  derivative  of  equation  (A3)  with 
respect  to  T  we  get  dP/dT  —  pk  since 

E/c°nf  f  Cconf(T).  (A4) 

Finally,  however  the  second  and  fourth  terms  of 
equation  (A2)  require  a  dF/dE  term.  From  the  simu¬ 
lation,  we  have  no  means  of  evaluating  this  term  and  a 
priori  knowledge  of  V  —  V(T)  is  unlikely.  A  functional 
fit  of  V  =  V(T)  could  be  determined  as  the  simulation 
proceeds  but  would  be  counter-productive  to  our  overall 
goal  of  determining  a  Hugoniot  state  from  a  single 
simulation  of  reasonable  length.  The  choice  of  x„  —  T 
in  equation  (2)  therefore  seems  less  practical. 
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~2dT 
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dE 


(C/conf) 


dE 


(EE) 
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Note  that  analogous  to  the  derivation  of  equation  (7), 
we  consider  the  quantities  in  this  expression  to  be 
instantaneous  values  that  depend  only  on  the  current 
configuration.  The  first  term  on  the  right-hand  side  of 
equation  (Al)  can  usually  be  analytically  determined 


References 

[1]  Touret,  J.,  and  van  den  Kerkhof,  A.  M.,  1986,  Physica 
B  &  C,  139,  834. 

[2]  Kortbeek,  P.  J.,  Seldam,  C.  A.  T.,  and  Schouten,  J.  A., 
1990,  Molec.  Phys.,  69,  1001. 

[3]  Duan,  Z.,  Moller,  N.,  and  Weare,  J.  H.,  1996,  Geochim. 
Cosmochim.  Acta,  60,  1209. 

[4]  Fickett,  W.,  and  Davis,  W.  C.,  1979,  Detonation 
(Berkeley:  University  of  California  Press). 

[5]  Ree,  F.  H„  1984,  J.  chem.  Phys.,  81,  1251. 

[6]  Fried,  L.  E.,  2001,  Cheetah  3.0  User’s  Manual,  (manuscript 
number  UCRL-MA-1 17541,  Revision  3)  (Livermore,  CA: 
Lawrence  Livermore  National  Laboratory). 

[7]  Swanson,  D.  R.,  Mintmire,  J.  W.,  Robertson,  D.  H., 
and  White,  C.  T.,  2000,  Chem.  Phys.  Reps.,  18,  1871. 


3322 


J.  K.  Brennan  and  B.  M.  Rice 


[8]  Robertson,  D.  H„  Brenner,  D.  W.,  and  White,  C.  T., 
1991,  Phys.  Rev.  Lett.,  67,  3132. 

[9]  Maillet,  J.  B.,  Mareschal,  M.,  Soulard,  L., 
Ravelo,  R.,  Lomdahl,  P.  S.,  Germann,  T.  C.,  and 
Holian,  B.  L.,  2000,  Phys.  Rev.  E,  63,  016121. 

[10]  Brennan,  J.  K.,  and  Rice,  B.  M.,  2002,  Phys.  Rev.  E,  66, 
021105. 

[11]  Shaw,  M.  S.,  2002,  12th  Detonation  Symposium,  San 
Diego,  11-16  August;  2001,  Shock  Compression  of 
Condensed  Matter,  AIP  Conference  Proceedings, 
Atlanta,  24-29  June,  edited  by  M.  D.  Furnish,  N.  N. 
Thadhani  and  Y.  Horie  (Melville,  NY:  AIP). 

[12]  Erpenbeck,  J.  J.,  1992,  Phys.  Rev.  A,  46,  6406. 

[13]  Rice,  B.  M.,  Mattson,  W.,  Grosh,  J.,  and  Trevino, 
S.  F„  1996,  Phys.  Rev.  E,  53,  611. 

[14]  Karo,  A.  M.,  Hardy,  J.  R.,  and  Walker,  F.  E.,  1978, 
Acta  Astronaut.,  5,  1041. 

[15]  Powell,  J.  D.,  and  Batteh,  J.  H.,  1978,  J.  appl.  Phvs.,  49, 
3933. 

[16]  Powell,  J.  D.,  and  Batteh,  J.  H.,  1979,  Phvs.  Rev.  B,  20, 
1398. 

[17]  Powell,  J.  D.,  and  Batteh,  J.  H.,  1980,  J.  appl.  Phvs.,  51, 
2050. 

[18]  Holihan,  B.  L.,  Hoover,  W.  G.,  Moran,  B.,  and 
Straub,  G.  K„  1980,  Phys.  Rev.  A,  22,  2798. 

[19]  Tsai,  D.  H.,  and  Trevino,  S.  F.,  1984,  J.  chem.  Phys.,  81, 
5636. 

[20]  Tsai,  D.  H.,  1990,  Chemistry  and  Physics  of  Energetic 
Materials,  edited  by  S.  N.  Bulusu  (Dordrecht:  Kluwer). 

[21]  Brenner,  D.  W.,  1992,  Shock  Compression  of  Condensed 
Matter  (Amsterdam:  Elsevier). 

[22]  Robertson,  D.  H.,  Brenner,  D.  W.,  Elert,  M.  L.,  and 
White,  C.  T.,  1992,  Shock  Compression  of  Condensed 
Matter  (Amsterdam:  Elsevier). 

[23]  White,  C.  T.,  Robertson,  D.  H.,  Elert,  M.  L., 
and  Brenner,  D.  W.,  1992,  Microscopic  Simulations 
of  Complex  Hydrodynamic  Phenomena  (New  York: 
Plenum). 

[24]  Brenner,  D.  W.,  Robertson,  D.  H„  Elert,  M.  L.,  and 
White,  C.  T„  1993,  Phys.  Rev.  Lett.,  70,  2174. 

[25]  White,  C.  T.,  Sinnott,  S.  B„  Mintmire,  J.  W., 
Brenner,  D.  W.,  and  Robertson,  D.  H.,  1994,  Int.  J. 
quantum  Chem.  Symp.,  28,  129. 

[26]  Rice,  B.  M.,  Mattson,  W.,  Grosh,  J.,  and  Trevino, 
S.  F.,  1996,  Phys.  Rev.  E,  53,  623. 

[27]  Alper,  H.  E.,  Abu-Awwad,  F.,  and  Politzer,  P.,  1999, 
J.  phys.  Chem.  B,  103,  9738. 

[28]  Bedrov,  D.,  Smith,  G.  D.,  and  Sewell,  T.  D.,  2000, 
J.  chem.  Phvs.,  112,  7203. 

[29]  Bedrov,  D.,  Smith,  G.  D.,  and  Sewell,  T.  D.,  2000, 
Chem.  Phys.  Lett.,  324,  64. 

[30]  Bedrov,  D.,  Ayyagari,  C.,  Smith,  G.  D.,  Sewell,  T.  D., 
Menifoff,  R.,  and  Zaug,  J.  M.,  2001,  J.  comput -aided 
mat.  Design,  8,  77. 

[31]  Bunte,  S.  W.,  and  Sun,  H.,  2000,  J.  phvs.  Chem.  B,  104, 
2477. 


[32]  Bunte,  S.  W.,  and  Miller,  M.  S.,  2001,  Army  Research 
Laboratory  Technical  Report,  ARL-TR-2496  (Aberdeen 
Proving  Ground,  MD:  US  Army). 

[33]  Smith,  G.  D.,  and  Bharadwaj,  R.  K.,  1999,  J.  phvs. 
Chem.  B,  103,  3570. 

[34]  Sorescu,  D.  C.,  Rice,  B.  M.,  and  Thompson,  D.  L., 

1997,  J.  phys.  Chem.  B,  101,  798. 

[35]  Sorescu,  D.  C.,  Rice,  B.  M.,  and  Thompson,  D.  L., 

1998,  J.  phys.  Chem.  A,  102,  8386. 

[36]  Sorescu,  D.  C.,  Rice,  B.  M.,  and  Thompson,  D.  L., 

1998,  J.  phys.  Chem.  A,  102,  8386. 

[37]  Sorescu,  D.  C.,  Rice,  B.  M.,  and  Thompson,  D.  L„ 

1999,  J.  phys.  Chem.  A,  103,  989. 

[38]  Sorescu,  D.  C.,  and  Thompson,  D.  L.,  1999,  J.  phvs. 
Chem.  B,  103,  6774. 

[39]  Sorescu,  D.  C.,  Rice,  B.  M.,  and  Thompson,  D.  L., 

1999,  J.  phys.  Chem.  B,  103,  6783 

[40]  Sorescu,  D.  C.,  Rice,  B.  M.,  and  Thompson,  D.  L., 

2000,  J.  phys.  Chem.  B,  104,  8406. 

[41]  Sorescu,  D.  C.,  Boatz,  J.  A.,  and  Thompson,  D.  L„ 

2001,  J.  phys.  Chem.  A,  105,  5010. 

[42]  Weisstein,  E.  W.,  Eric  Weisstein’s  World  of  Mathematics 
(http://mathworld.wolfram.com/ Root-Finding  Algorithm, 
html). 

[43]  Kofke,  D.  A.,  1993,  Molec.  Phys.,  78,  1331. 

[44]  Kofke,  D.  A.,  1993,  J.  chem.  Phys.,  98,  4149. 

[45]  Mehta,  M„  and  Kofke,  D.  A.,  1994,  Chem.  Eng.  Sci., 
49,  2633. 

[46]  Allen,  M.  P„  and  Tildesley,  D.  J.,  1987,  Computer 
Simulation  of  Liquids  (Oxford:  Oxford  University  Press). 

[47]  Fried,  L.  E.,  and  Howard,  W.  M.,  1998,  J.  chem.  Phvs., 
109,  7338. 

[48]  Frenkel,  D.,  and  Smit,  B.,  2002,  Understanding 
Molecular  Simulation  (San  Diego:  Academic  Press). 

[49]  Reed,  T.  M.,  and  Gubbins,  K.  E.,  1973,  Applied 
Statistical  Mechanics  (New  York:  McGraw-Hill). 

[50]  Lofthus,  A.,  and  Krupenie,  P.  H.,  1977,  J.  phys.  Chem. 
ref.  Data,  6,  113. 

[51]  Moore,  C.  E.,  and  Gallagher,  J.  W.,  1993,  Tables  of 
Spectra  of  Hydrogen,  Carbon,  Oxygen  Atoms  and  Ions, 
CRC  Series  in  Evaluated  Data  in  Atomic  Physics  (Boca 
Raton:  Chemical  Rubber). 

[52]  Melchionna,  S.,  Ciccotti,  G.,  and  Holian,  B.  L.,  1993, 
Molec.  Phvs.,  78,  533. 

[53]  Johnson,  J.  K.,  Panagiotopoulos,  A.  Z.,  and  Gubbins, 
K.  E.,  1994,  Molec.  Phvs.,  81,  717. 

[54]  Smith,  W.  R.,  and  Triska,  B.,  1994,  J.  chem.  Phvs.,  100, 
3019. 

[55]  Johnson,  J.  K.,  1999,  Adv.  chem.  Phys.,  105,  461. 

[56]  McQuarrie,  D.  A.,  1976,  Statistical  Mechanics 

(New  York:  Harper). 

[57]  Kristof,  T.,  and  Liszi,  J.,  1998,  Molec.  Phys.,  94,  519. 

[58]  Ray,  J.  R.,  and  Graben,  H.  W.,  1986,  Phys.  Rev.  A,  34, 
2517. 

[59]  Smith,  W.  R.,  and  Lisal,  M.  L.,  2002,  Phys.  Rev.  E,  66, 
011104. 


NO.  OF 

COPIES  ORGANIZATION 


1  DEFENSE  TECHNICAL 
(PDF  INFORMATION  CTR 
ONLY)  DTICOCA 

8725  JOHN  J  KINGMAN  RD 
STE  0944 

FORT  BELVOIR  VA  22060-6218 

1  US  ARMY  RSRCH  DEV  & 
ENGRG  CMD 
SYSTEMS  OF  SYSTEMS 
INTEGRATION 
AMSRD  SS  T 
6000  6TH  ST  STE  100 
FORT  BELVOIR  VA  22060-5608 

1  DIRECTOR 

US  ARMY  RESEARCH  LAB 
IMNE  ALC  IMS 
2800  POWDER  MILL  RD 
ADELPHI  MD  20783-1197 

3  DIRECTOR 

US  ARMY  RESEARCH  LAB 
AMSRD  ARL  Cl  OK  TL 
2800  POWDER  MILL  RD 
ADELPHI  MD  20783-1197 


ABERDEEN  PROVING  GROUND 


1  DIR  USARL 

AMSRD  ARL  Cl  OK  TP  (BLDG  4600) 


NO.  OF 

COPIES  ORGANIZATION 

ABERDEEN  PROVING  GROUND 


5  DIR  USARL 

AMSRD  ARL  WM  BD 
J  BRENNAN 


